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We present here a detailed numerical study of the dynamical behaviour of ‘soft’ uncompressed 
grains in a granular chain where the grains interact via the intrinsically nonlinear Hertz force. It 
is well known that such a chain supports the formation of solitary waves (SW). Here, however, 
the system response to the material properties of the grains and boundaries is further explored. In 
particular, we examine the details of the transition of the system from a SW phase to an equilibrium¬ 
like (or quasi-equilibrium) phase and for this reason we ignore the effects of dissipation in this study. 

We find that the soft walls slow the reflection of SWs at the boundaries of the system, which in 
turn slows the journey to quasi-equilibrium. Moreover, the increased grain-wall compression as the 
boundaries are softened results in fewer average grain-grain contacts at any given time in the quasi¬ 
equilibrium phase. These effects lead to increased kinetic energy fluctuations in the short term in 
softer systems. We conclude with a toy model which exploits the results of soft-wall systems. This 


toy model supports the formation of breather-like 
energy in desired places in the granular chain. 

PACS numbers: 45.70.-n, 05.45.Yv 

I. INTRODUCTION 

Following the pioneering works of Nesterenko [IHI3], 
a tremendous amount of experimental, analytical, and 
numerical studies into the dynamics of granular chains 
have emerged in recent years (see e.g. Refs. [MHH]). It 
is well established that, in general, any velocity pertur¬ 
bation to an end grain in a ID granular chain without 
any initial pre-compression, i.e. an unloaded chain, will 
propagate through the system as a nondispersive bundle 
of energy, or solitary wave (SW) [20]. This SW devel¬ 
ops within approximately ten grain diameters [22], and 
it travels through the chain with a constant speed which 
scales with SW amplitude, as is characteristic for nonlin¬ 
ear waves [20] . This SW also has a well-defined and fixed 
width that is sensitive to the grain shape m, and equal 
to approximately five grain diameters for grains with el¬ 
lipsoidal contact geometries such as spheres mm- 
The SW arises from the contact potential between 
grains described by the Hertz law [24], V ^ , where 

(5 > 0 is the grain-grain overlap and n = 5/2 for spher¬ 
ical grains. The intrinsically nonlinear Hertz potential 
possesses no quadratic term, and thus one cannot have 
harmonic motion in the unloaded granular chain. This 
manifests itself in the absence of acoustic propagation 
in the system [25], a fact referred to by Nesterenko and 
others as the “sonic vacuum” [T]. Adding compression 
into the chain effectively introduces a harmonic term to 
the potential, and the effects of this term on the exis- 
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entities and may therefore be useful for localizing 


tence and stability of SWs has been extensively studied 
elsewhere, see e.g. Refs. and [2^ 

The SW supported by the granular chain is unlike soli- 
tons found in other (i.e. continuum) physical systems: it 
is a unique SW specific to discrete lattices with purely 
nonlinear interactions. The unloaded chain’s discrete na¬ 
ture has the possibility of broken grain-grain contacts, 
and this crucial feature means that, in general, SWs are 
not perfectly preserved. Collisions of SWs with other 
SWs has been extensively probed in the literature, see 
e.g. Refs. [27l [28] and [29l SWs in the granular chain are 
also not preserved in collisions with boundaries [11130]- 
[32]. Rather, the SW breaks up and reforms in the col¬ 
lision process [28], resulting in the partial destruction of 
the original SW and the birth of secondary solitary waves 
(SSWs) which are much smaller in magnitude, but pos¬ 
sess the same spatial extent as the parent SW [33] . 

SSWs were first reported in numerical studies in 
Refs. |3l and |34| and later confirmed experimentally by 
Job et al. m- It is now known that the presence of 
boundaries in these systems generally leads to the forma¬ 
tion of SSWs with various energies which end up moving 
at a distribution of velocities related to the energy car¬ 
ried in the SSW. In the absence of dissipation and at 
sufficiently long time after an initial energy perturbation 
to the system, the rates of breakdown and creation pro¬ 
cesses of SSWs balance and the chain reaches a steady 
state which has been referred to in the literature as the 
quasi-equilibrium (QEQ) phase [29ll32] . Here we inves¬ 
tigate how the system’s journey to QEQ can be tuned 
by varying the material parameters of the granular chain 
system, and the effects of introducing an inertial mis¬ 
match at the boundary on the onset and rate of relax¬ 
ation to the QEQ phase. We will ignore the effects of dis- 



2 


sipation that are invariably present in a system made up 
of macroscopic elastic objects. A dissipation free system 
should be regarded as a first and necessary step in reach¬ 
ing driven-dissipative systems with QEQ-like behavior. 

Here we extend the work of Ref. [35] by considering the 
effects of either softening the boundary walls or the whole 
system on the dynamics of the unloaded granular chain 
held between fixed symmetric boundaries. We are par¬ 
ticularly concerned with the effects of soft boundaries on 
the rate of relaxation as the system slips into QEQ. As 
we show below, the time delay introduced by the soft 
walls in the rebounding of SWs at the boundaries in¬ 
creases with the wall softness. This time delay has two 
main effects on the system dynamics: in the short term, 
soft-wall systems exhibit larger kinetic energy fluctua¬ 
tions than corresponding systems with homogeneous or 
hard walls; in the long term, the transition to the QEQ 
phase is slowed down by incorporating soft boundaries 
into the system. We suggest that, in soft-wall systems, 
the journey from the initial SW state to the QEQ phase 
can be understood in terms of a transition from a SW 
state to an ergodic phase where the ensemble average of 
grain energy is equal to the time average. 

One key interest in the QEQ phase is with possible en¬ 
ergy harvesting and localization applications. Recent in¬ 
terest in these systems has probed metastable special re¬ 
gions such as so-called discrete breathers [361438] and hot 
and cold spots that are sustained for extended times [29] , 
where the kinetic energy can be localized along well- 
defined locations. Here we use insights from results of 
soft-wall systems to discuss a toy model which may be 
useful for locally trapping energy. 


II. METHODS 

We consider N monodisperse spherical grains in ID, 
held between two fixed symmetric walls, under zero pre¬ 
compression. Initially the grains are barely touching, but 
as any two grains in the chain make contact they inter¬ 
act via a contact potential given by the Hertz law [24] . 
V{Si^i+i) = where Si^i+i = {Ri+i - Ri) - 

— Xi) is the grain-grain overlap with Xi denoting 
the absolute position of grain i with respect to the ori¬ 
gin, and Ri its radius. The proportionality constant is 
given by 


^ 2,2 + 1 


2 / RjRj^i 

^Di^i^i \J Ri Ri-\-i ’ 


( 1 ) 


where the prefactor is related to the Young’s mod¬ 

ulus Y and Poisson’s ratio a of the grains by 
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We consider here chains in which all the grains are com¬ 
prised of the same material, and thus = (7^+1 and 


Yi = Vi. It is then possible to express the equa¬ 
tions of motion for any grain (except the end grains) in 
the chain as follows: 
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mxi = - 


. 3/2 . 3/2 


(3) 


where m is the grain mass, Xi denotes the second deriva¬ 
tive of the absolute position with respect to time, and 
i = 2... Y — 1. To implement the fixed walls, one takes 
Rj ^ oc as the radius of the wall, ensuring that the 
boundary is much larger than any grain in the chain and 
will therefore not move during the simulation, while si¬ 
multaneously relaxing the condition that the wall must 
be flat. In this way, the potential that develops at the 
interface between either the first or the last grain in the 


chain and the wall is given by: V = {2/bDi^j)\/R5^JJ‘^ 
where R is the grain radius for the end grains. 

In this work we refer to a chain in which the grains 
and walls are comprised of the same material as a ho¬ 
mogeneous chain. One requires only a single force pref¬ 
actor to fully characterize the system in this case, since 
Di,i+i = Dij = D. The potential that develops between 
an end grain and wall is then a factor of \/2 larger than 
the grain-grain potential in the homogeneous systems. 
The chains in which the grain material differs from the 
wall material are referred to in this work as inhomoge¬ 
neous systems. In this case, one requires two force pref¬ 
actors to characterize the system: one for the grain-grain 
interactions, = Dg^ and one for the grain-wall in¬ 

teractions Dij = D^. Clearly, Dg = in homogeneous 
granular systems. 

We use the velocity-verlet algorithm [39] with a time 
step of dt = 10“^ /is to integrate the equations of motion. 
Initially all of the grains are at rest, and the left end grain 
is given a velocity perturbation directed into the chain at 


t = 0 /is of magnitude vq = 9.899 x 10 


i//is in all 


cases. The system is then left to evolve in time according 
to the equations of motion. Typically, chains were com¬ 
prised of Y = 38 spherical grains, and all grains were 
taken to be R = 6 mm in radius. The choice of system 
size was linked to probing chains which may be experi¬ 
mentally accessible, and we show in Sec. HI A that this 
is a suitable choice. We considered a variety of materials 
for the homogeneous systems whose material parameters 
can be found in the appendix, sorted according to the 
product of the mass density p and prefactor D (vide in¬ 
fra). These materials range from as hard as diamond 
{pD = 0.00415 (mm/yas)“^) to as soft as silicone rubber 
{pD = 52.4 (mm//is)“^). 

Eor inhomogeneous systems, all chains were taken to be 
comprised of grains with Y = 193 GPa, a = 0.305, and 
p = 7.82 mg/mm^, consistent with material properties 
of stainless steel. This corresponds to a grain-grain force 
prefactor of Dg = 0.0071 (mm/is^/mg). Keeping the 
grain material fixed, we then symmetrically vary the wall 
material parameters to elucidate the effects of softening 
the boundaries on the system dynamics. 

Eor the inhomogeneous systems, the values were 
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chosen to vary over a suitable range, thus some are con¬ 
trived and do not correspond to any particular material. 
We considered = 0.0041 (diamond walls), 0.0071 
(stainless steel walls), 0.05, 0.5, 5, 10, 19 (~ rubber 
walls), 30, 50, and 75 (mm/is^/mg). A peculiar fea¬ 
ture of these systems is that the notion of “hard” and 
“soft” walls is not absolute, but is relative to the mate¬ 
rial comprising the grains in the chain, and the walls can 
be considered soft only if > Dg. For example, the 
stainless steel systems considered here begin to exhibit 
soft-wall behaviour only once > 0.0071. 


III. NUMERICAL RESULTS 

After initiating a delta-function velocity perturbation 
to the left end grain at t = 0 /is, it takes several grains for 
the SW to develop in the system, as illustrated by kinetic 
energy density plots in Figs. and A small amount 
of kinetic energy remains at the left wall due to the left 
end grain rebounding off its neighbour to transfer the 
kinetic energy of the initial solitary wave (ISW), down 
the chain. Likewise, after the ISW reaches and emerges 
from the right wall, secondary solitary waves (SSWs) are 
formed at the right wall from a similar rebounding of the 
right end grains. These energy processes occur in both 
the inhomogeneous and homogeneous systems and are 
visible for a variety of cases in the early time dynamics 
of Figs. and 

The most striking feature of Fig. is the emergence 
of a time delay in the ISW rebounding at the boundary 
in the inhomogeneous systems as one progressively soft¬ 
ens the walls. This appears as a horizontal line segment 
near grains 1 and 38, and is absent in the hard-wall and 
homogeneous systems. In contrast to homogeneous sys¬ 
tems, soft wall systems take a notable amount of time 
1600 yus for = 19) to rebound from the compres¬ 
sion of the end grain as the ISW collides with the far wall. 
This slow time scale associated with interactions at the 
walls and a much faster time scale associated with SW 
interactions within the chain (i.e. away from the wall) 
leads to an inertial mismatch between the grains and the 
walls, or the appearance of “sticky” walls in the inhomo¬ 
geneous systems. One finds that the softer the wall is 
made to be, the longer the time it takes to re-emit the 
SW from the boundary. 

According to the virial theorem, a system with poten¬ 
tial energy V ^ 5^ will partition the total system energy 
as n/(n + 2) kinetic energy and 2/(n + 2) potential en¬ 
ergy. For n = 5/2, this leads to 5/9 of the total energy 
as kinetic energy, and the remaining 4/9 as potential en¬ 
ergy. Most of the system’s energy is held as potential 
energy while the end grains are compressed against the 
soft walls for extended periods of time. Thus the soft- 
wall systems make larger SSWs to compensate, in an 
attempt to partition the system’s total energy as appro¬ 
priate. This process is particularly visible by comparing 
kinetic energy density plots in Fig. for = 0.007 


(homogeneous) and = 19: the SSW emitted at the 
right wall near 0.0025 s in the inhomogeneous soft-wall 
system has a much larger energy content, thus appearing 
“warmer” than that emitted in the homogeneous case. 

Another interesting feature of Fig. are the “bending 
waves” which appear as curved line segments (e.g. be¬ 
ginning at about 0.004 s for = 19) and occur more 
prominently in the soft-wall systems than homogeneous 
systems. Any time grains break contact, a gap is formed 
in the system causing the SW to lose its coherence and 
partially break. Thus some fraction of the SW kinetic 
energy located in broken-contact grains becomes decou¬ 
pled from the main part of the wave, forming new small 
SSWs. Propagation of a SSW is slowed by gaps between 
grains; until the grain with peak kinetic energy catches 
up with the rest of the chain, the progress of the SSW 
cannot continue. As inter-grain gaps are closed, the SSW 
picks up speed until reaching its natural velocity, which 
is determined by the Hertz potential exponent, wave am¬ 
plitude, and grain material properties [20]. Thus many 
SSWs appear to ‘curve’ in the kinetic energy density plot. 
We speculate that soft-wall systems facilitate the opening 
of larger gaps in the system, leading to more prominent 
bending of the SSWs. 

An increased tendency to open larger and more fre¬ 
quent gaps in the system intuitively results from the soft 
walls being more easily compressed by the harder grains. 
In fact, the first and last grains always remain in con¬ 
tact with the soft walls during the rebounding of SWs 
and the emission of SSWs at the boundaries. Once the 
last grain initially makes contact with the right wall in 
a soft-wall system, it never breaks contact for the dura¬ 
tion of the simulation. This is also true for the left end 
grain, after rebounding off the second grain. The sus¬ 
tained grain-wall compression facilitates the opening of 
gaps near the ends of the chain, which are evident from 
an analysis of individual grain kinetic energies over time. 
For the steel grain/rubber wall = 19) case for ex¬ 
ample, it takes a significant amount of time for grain 38 
to pass its kinetic energy to grain 37, and also for 37 
to pass to 36, etc. as the ISW is rebounding from the 
right wall. This is because the end grains have become 
disconnected from each other due to the opening of a no¬ 
ticeably large gap at the right end of the system as grain 
38 remains compressed into the right wall. Consequently, 
end grains travel with a constant velocity before striking 
their neighbours: grain 38 travels ~ 4 /rm before hitting 
grain 37, and grain 37 travels 10 times this distance 
before colliding with grain 36 in this system. 


A. Transition to quasi-equilibrium 

The build up of SSWs from inter-grain gaps and ex¬ 
tended wall-grain contacts eventually destroys the phase 
dominated by the ISW, and the system transitions into a 
new phase with many small SSWs breaking and reform¬ 
ing very quickly. This can be seen in the emergence of 
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FIG. 1. (Color online) Kinetic energy density plots for the N = 38 inhomogeneous stainless steel granular chains with varying 
wall softness D^. For each panel (a) - (d), Dg = 0.0071. The kinetic energy values are normalized by the total system energy 
and presented using a logarithmic scale. The SW speed is set by the grain material, see Fig.[^ and is the same for all systems 
presented here. Notice the progressive delay that develops in the rebounding of SWs at the boundaries as they are made 
increasingly softer. 


noisy energy fluctuations in Fig. at longer times. 

This g'^asz-equilibrium (QEQ) phase for the initial con¬ 
ditions used here is marked by a Maxwell-Boltzmann dis¬ 
tribution of particle velocities, but where the equipar- 
titioning of energy among all grains is absent. Since 
the ISW is partially destroyed in collisions with system 
boundaries- which, as a byproduct slowly produce SSWs 
in the vicinity of the collision region- the energy content 
of the ISW decreases over time, and energy is partially 
distributed among other grains in these unbiased sys¬ 
tems. The motion of individual grains in one-dimensional 
systems is highly correlated, hence groups of particles of¬ 
ten oscillate together. The fastest grain in the chain at 
any instant will either be part of the ISW, in which case 
it will have a large share of the system’s kinetic energy, 
or it will have only a small share of the kinetic energy of 
several SSWs in the QEQ phase. Therefore to visualize 
the transition from the unstable phase made of various 
SWs and SW-like structures to the QEQ phase, one can 
consider simply the speed of the fastest grain in the chain. 

We present a series of plots in Eig. which depict the 
absolute value of the velocity of the fastest moving grain 
in the chain, \vmax\^ as a function of time for several sys¬ 
tems. Initially, \vmax\ oscillates about a maximal value 
as the ISW travels through the granular chain. Peri¬ 


odically, when the ISW reaches a boundary most of the 
system’s energy converts to stored potential energy in the 
walls, and \vmax\ drops to a small value. This happens 
only briefly in the homogeneous systems since the time it 
takes for grains to transfer kinetic energy to each other is 
approximately equal to the duration of rebounding of the 
ISW off the wall. However, there is an additional time 
delay in the rebounding of SWs at the soft wall in the 
inhomogeneous systems, and the interplay between these 
time scales leads to the emergence of a beating pattern 
in the soft-wall systems, apparent in Eigs.j^a) and|^b). 
In this case, the inherent time scale set by the material 
properties of the walls is not equivalent to the grain-grain 
energy transfer time, as it is in the homogeneous systems. 
It is possible that one might be able to exploit this fre¬ 
quency mismatch and beating response by constructing 
piezoelectric systems with a targeted frequency response. 


Some time after the initial velocity perturbation, 
I'^maxl relaxes to noisy oscillations about a much smaller 
value, denoting the onset of the QEQ phase. One no¬ 
ticeable feature of Eigs. ia) and l^b) is the delay of the 
relaxation process introduced by softening the walls. We 
model this phase transition with a sigmoidal error func- 
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(a) D = 0.001 (diamond) 



time (s) 



(d)D = 38 (rubber) 



time (s) 


FIG. 2. (Color online) Kinetic energy density plots for the homogeneous N = 38 granular chains comprised of varying material 
parameters. Kinetic energy values were normalized by the total system energy and presented on a logarithmic scale. The 
time scale has been adjusted to capture the first four rebounds of the ISW at the right wall. Recall for homogeneous systems 
Dw = Dg = D. The SW possesses a speed which is set by the softness of the grains, and progressively softer systems {pD ^ 1) 
correspond to smaller wave speeds. 


tion of the form: 


In this way one can extract the relaxation time cr[/is], as 
well as the midpoint of relaxation to QEQ r* [//s]. A com¬ 
parison of (j and are presented in Figs.Qa) andj^b), 
respectively, for inhomogeneous systems, and Figs.^j^c) 
and l^d), respectively, for homogeneous systems. This 
two-state transition model makes sense for the homo¬ 
geneous cases, where \vmax\ exhibits a clear changeover 
from SW state to the QEQ state, e.g. Figs. 3(c) and 
3(d). 

However, it is more difficult to implement for the soft- 
wall systems, where energy is locked-up in the walls as 
potential energy for extended periods of time, e.g. Figs. 
3(a) and 3(b). These gaps in the trend of \vmax \ prohibit 
robust fitting of Eq. and lead to poor goodness of fits 
^ 0.5 — 0.6). Unless a better statistical measure of 
the transition can be found, the exact values of a and 
for the inhomogeneous case must be considered impre¬ 
cise, although the overall trend is qualitatively correct. 

Results for the inhomogeneous stainless steel systems, 
Figs.j^a) and|^b), illustrate that one can not only slow 


the rate of relaxation of the system to QEQ beyond its 
corresponding homogeneous value, but one can also phys¬ 
ically delay the onset of the journey from the SW state 
to the QEQ state by softening the walls of the system. 
We reiterate here that “soft” is a relative term in these 
systems, and one can define the walls to be soft only if 
g • T^he important point is that the v^alls can be 
made arbitrarily soft, and thus it should be possible to 
somewhat tune or delay the transition to QEQ as needed 
in a real physical system. 

For the homogeneous systems, Figs.|^c) andj^d), we 
find a strict power-law dependence of both a and r* upon 
the parameter pD^ where p is the grain/wall density and 
D = Dg = is the force prefactor, introduced in Eq. 

In particular, one finds that as the grains and walls are 
simultaneously softened in homogeneous systems, both 
the relaxation time and onset of relaxation are shifted to 
higher values. 

To further characterize the granular chains, we con¬ 
sidered the fluctuations in kinetic energy from the value 
predicted by the virial theorem. In particular, the fluc¬ 
tuation in the system’s total kinetic energy (Tk was com¬ 
puted as the relative root mean squared deviation over 


= vx ( erfc 
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(b) Dw=50 



— y = a.erfc(b(x-c))+d 
a =1.68x10'^ 
b = 9.59x10'^ 
c = 1.48x10^ 
d = 2.96x10'^ 

R^ = 0.6837 


10000 20000 30000 40000 50000 60000 
time (jLis) 


(d) D = 0.02 (pyrex) 

—I - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 

— y = a.erfc(b(x-c))+d 



10000 20000 30000 40000 50000 
time (ps) 


FIG. 3. (Color online) Speed of the fastest moving grain in the N = 3S granular chain as a function of time, (a) and (b) 
correspond to inhomogeneous stainless steel (Dg = 0.007) granular chains with varying wall softness Dw (c) and (d) correspond 
to homogeneous {D^ = Dg = D) granular chain systems. 


intervals m which run from time tm ^ tm t: 

In the above equation, r is the considered time inter¬ 
val (100 jiis here), corresponds to the number of data 
points in this interval, and the mean of the system’s ki¬ 
netic energy (K)^ is a constant equal to the value pre¬ 
dicted by the virial theorem, and not the actual average 
(K) of the interval. Since all systems were given the same 
initial velocity perturbation, the total energy of homoge¬ 
neous systems depends upon the material properties of 
the granular chain itself, and {K)^ is therefore different 
for each of the homogeneous cases. 

We compute the average fluctuation in the system’s 
kinetic energy as the average of the intervals m' from the 
start of the simulation, 

1 ^ 

«(iM) = 'D. (6) 

m' = l 

These results are shown in Fig. |^a) for stainless steel 
chains with varying wall softness, as well as for an ar¬ 


ray of homogeneous systems in Fig. ^c). For the in¬ 
homogeneous systems, we see that ^{t) exhibits a series 
of peaks which are modulated by an overall exponential 
decay at early times, and at much longer times, K{t) re¬ 
laxes to a long tail with algebraic dependence (we find, 
e.g., oo) ^ 1/t for the N = 3S stainless steel chain 

with = 5). Deep into QEQ, the average fluctuation in 
the system’s total kinetic energy levels off at a constant 
value which is insensitive to the wall softness parameter 
Dw, and for chains of length = 38 is ^ 15% of the 
kinetic energy value predicted by the virial theorem. 

We performed additional simulations of chains of vary¬ 
ing length N for homogeneous stainless steel systems, as 
well as inhomogeneous chains with stainless steel grains 
and rubber walls {D^ = 19). The long-term average fluc¬ 
tuation is sensitive to the system size, K{t oc) ^ 
(approximately) for N ^ 1 in both cases, as shown in 
Figs.[^b) and 5(d). Thus kinetic energy partitions more 
evenly for large systems. 

B. Long-term behaviour 

Deep in the QEQ phase, for this study taken to be the 
last 10% of the simulation, an analysis of the grain veloc- 
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logio[pD^ (mm/gs)-^] 



FIG. 4. (Color online) Results of complementary error function fits to \vmax\{t) from Fig. with logarithmic scales on both 
axes. Plots (a) and (b) display, respectively, the dependence of the relaxation time a, and the midpoint of relaxation to QEQ 
T*, on the wall material parameter pD^. (c) and (d) illustrate how the relaxation time and midpoint of relaxation, respectively, 
depend on the grain material parameter pD in homogeneous systems. Solid lines in (c) and (d) illustrate fits to the homogeneous 
results. 


ities for both homogeneous and inhomogeneous systems 
confirms a Maxwell-Boltzmann distribution, as shown in 
Fig-i substantiating our choice of system size N was 
not too small. Analysis of the temperature T obtained 
via Maxwell-Boltzmann fits of the velocity distributions 
for each system reveals that T is approximately con¬ 
stant with respect to the wall softness parameter 
as implied by Fig. ib). Moreover, using the temper¬ 
ature from all the stainless steel and soft-wall systems, 
the total kinetic energy as determined by NksT{2^ is in 
agreement with the calculated average total kinetic en¬ 
ergy determined by the average of the grain velocities, 
e.g. {J2mvf/2). 

Further into the QEQ phase, the kinetic energy is fairly 
evenly distributed among the grains in our systems, thus 
any evidence of “cold” or “hot” spots (such as e.g. in 
the vicinity of soft walls) have been erased at this point. 
It turns out that while equipartitioning of energy is ex¬ 
pected to be absent in QEQ, it is more difficult to quan¬ 
tify the violation of the equipartition theorem than one 


might imagine [20] . The presence of hot and cold spots 
from figures similar to Figs. and [^have been observed 
in numerical simulations of granular systems in the QEQ 
phase m and is one of the characteristics which makes 
this state different from a conventional equilibrium state. 
We additionally note that the ergodicity of the system 
deep into QEQ was checked as in Ref. [35|, and both in¬ 
homogeneous and homogeneous systems were found to 
be ergodic. Soft walls have the global effect of reducing 
the average total number of grain-grain contacts in the 
long-term dynamics of the system. This is illustrated in 
Fig- [^a), which displays the average number of grain- 
grain and grain-wall contacts deep into QEQ for N = 3S 
stainless steel granular chains with varying wall softness 
parameter D^. Note that the overall magnitude of the er¬ 
ror bars remains relatively constant in Fig.j^a), thus the 
relative fluctuation in the number of contacts increases 
as a function of D^. 

A similar plot is shown in Fig. [^c) for homogeneous 
chains, and it is evident that the average total number 








FIG. 5. (Color online) (a) Average fluctuation in total kinetic energy for iV = 38 stainless steel granular chains with varying 
wall softness parameter Dw (b) Average long-term fluctuation oo) as a function of chain length N in the stainless steel 

granular chain with rubber walls (Dw = 19). (c) Average long-term kinetic energy fluctuation for the N = 38 homogeneous 
chains of varying material parameters, (d) Average long-term fluctuation as a function of system size for homogeneous stainless 
steel systems. Lines in (b) and (d) have slope of —1/2, indicating the approximate dependence K{t oc) ~ N~^'^. 


of contacts is independent of the material comprising the 
granular chain. Thus a softening of the boundaries in 
comparison to grains is required to reduce the average 
number of contacts. The fact that the end grains spend 
the whole simulation in close contact with the softer walls 
opens a lot of space between the walls, which is then 
distributed between the grains. Paradoxically, the much 
reduced number of average grain-grain contacts with soft 
walls does not speed the transition to QEQ, c.f. Figs. Eta) 
and [^b). This leads one to speculate that while the 
number of gaps increases, they are more transitory but 
occur more frequently. 

The unbroken contact between the end grains and walls 
gives rise to a steady force between them, even as they 
continue to transfer kinetic energy to neighbouring grains 
during the reflection of SWs and emission of SSWs at 
the boundaries. Despite this sustained force between 


the end grains and walls, the magnitudes of the forces 
are increasingly smaller in softer systems. Thus the in¬ 
terplay between increased duration of applied force and 
decreased magnitude of force leads to an overall smaller 
time-integrated force (or impulse) in the long term. Over 
a long time, one sees a linear dependence of the cumu¬ 
lative impulse, I{t) = , between the wall and 

end grain interfaces with time, as seen in the insets to 
Figs, [^b) and 7(d). 

It is evident from Fig. [^b) that for a given granular 
chain, the slope of the impulse curve decreases exponen¬ 
tially as the walls are made softer. Similarly, the slope of 
the impulse curve is seen to decrease in homogeneous sys¬ 
tems as the entire system is softened (i.e. as D increases), 
see Fig. IZtd)- (Note the — log scale of Fig. [^d).) Inspec¬ 
tion of this plot reveals there are anomalies to this rule, 
suggesting that the density of the grains also plays a role 
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FIG. 6. (Color online) (a) Maxwell-Boltzman fits to the normalized velocity distributions in the inhomogeneous stainless 
steel systems with varying wall softness parameter D^. (b) Comparison of average kinetic energy deep into quasi-equilibrium 
for stainless steel chains of varying wall softness. Circles are the total kinetic energy K = NksT obtained from the Maxwell- 
Boltzmann temperature, and squares are K — mv^ j2) from the equations of motion, (c) and (d) illustrate the corresponding 
plots for homogeneous systems. Note that a logarithmic scale has been used in (d). The average kinetic energy is necessarily 
different for each system since all systems were given the same initial velocity (not energy) perturbation. Thus, we expect to 
see an agreement between the two curves in (d), which overlap almost perfectly. 


in fully characterizing the “softness” of homogeneous sys¬ 
tems. Nevertheless, comparisons of Figs, [^b) and 7(d) 
reveal that homogeneous systems = Dg = D) have 
a much smaller impulse slope at the end grain/wall in¬ 
terface than an inhomogeneous soft-wall system with the 
same value. 

A linear dependence of the cumulative impulse on time 
implies that F(t) is a constant of motion in the QEQ 
phase. Thus the impulse slope is the average force, and 
the slopes of linear fits to the impulse agree well with 
the average forces calculated from simulation data deep 
into QEQ (i.e. from averaging late-time data obtained 
from directly integrating the equations of motion). Since 
the force of compression between any two grains is re¬ 
lated to the grain-grain overlap via the Hertz potential 
as F = (5/2)0(53/2, we can calculate the average over¬ 


lap between the end grain and adjacent wall through the 
slope of the impulse curve. Using the expressions from 
Sec. [TO and taking Rj ^ oc as the radius of the wall, 
it follows that the average overlap parameter at the end 
grain/wall interface ((5^/^) = I^wEavg/ a/ 6 where the grain 
radius is taken to be i? = 6 mm. 

We show the average end grain/wall overlap param¬ 
eter for homogeneous and inhomogeneous systems as a 
function of pD^ in Eig. Here we include the grain 
density for all cases to facilitate full characterization of 
the softness of the homogeneous systems. The data pre¬ 
sented in Eig. were calculated from the slopes of linear 
fits to the impulse curves shown in Eigs. [^b) and[^d). 
On the log-log plot, see Eig.[^ the homogeneous systems 
exhibit a linear relationship with the parameters shown; 
however, the inhomogeneous systems are well described 
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FIG. 7. (Color online) (a) Average number of grain-grain and grain-wall contacts deep into QEQ for the N = 38 stainless steel 
granular chain with varying wall softness, Dw The grains are no longer in contact when (5 < 0. (b) The slope of the cumulative 
impulse with respect to simulation time, as a function of wall softness. The impulse is perfectly linear in time, shown in the 
inset for large and small examples of Dw (c) Average number of grain-grain and grain-wall contacts deep into QEQ for the 
homogeneous N = 38 granular chains of varying hardness D. (d) The slope of the cumulative impulse with respect to D. Note 
the slope is reported as — log scale, since the very rapid drop to small values as D > 1. The impulse is again perfectly linear 
in time, shown in the inset for large and small examples of D. 


by a quadratic curve in the figure. We find that deep 
into the QEQ phase, the system attempts to distribute 
forces roughly equally among all interfaces in the system. 
Hence the overlap at the grain/wall interface is larger if 
one introduces soft walls, which allows this interface to 
achieve the same force of compression that exists between 
grain-grain interfaces. 


IV. A TOY MODEL 

Finally, we discuss a toy model granular system de¬ 
signed for locally trapping energy, and the results of some 
simulations of this system. In particular, we draw on the 
results of the previous section, exploiting the “sticky” en¬ 
ergy propagation associated with the large time scales set 
by soft materials. We consider in this section a granular 


chain with soft boundaries and also central soft-grain im¬ 
purities. The soft boundaries may not be necessary for 
facilitating the localization of energy at the center of the 
chain, but are implemented to ensure the system is made 
as soft as possible. Here we present the results of sim¬ 
ulations on soft-wall systems comprised of three or four 
soft grains located centrally within a chain of hard grains 
which have double the radius of the impurities, as shown 
schematically in Figs, [^a) andj^c). Changing the ra¬ 
dius of the impurity does not seem to have a significant 
effect on the results, though further analysis is required 
to fully characterize the role of the size of the impurity 
on the energy localization process. The main result is 
that there is a remarkable slowing down in the system in 
the vicinity of the soft impurities. The resulting energy 
objects are not necessarily so-called discrete breathers, 
which are generally described as time-periodic, localized. 
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logjQ[pD^ (mm/|Lis)'^] 

FIG. 8. (Color online) The average overlap parameter 
as a function of material softness pHw, where p is the grain 
density and Hw the force prefactor for the end-grain/wall in¬ 
terface. The average end-grain/wall overlap parameter was 
calculated from linear fits to the cumulative impulse. The 
intersection point marks the homogeneous stainless steel sys¬ 
tem. Note that the data is presented on a logarithmic scale 
for both axes, and the y-axis denotes the negative logarithm 
of Thus a smaller value on the y-axis corresponds to 

a larger overlap. The single square to the left of the inter¬ 
section point corresponds to the stainless steel system with 
harder diamond walls. 


and stable excitations in spatially extended discrete sys¬ 
tems. However, the localization of the energy on the 
impurity behaves effectively like a short-lived breather- 
breather, exhibiting an interesting natural frequency re¬ 
sponse resulting from the mutual repulsion between the 
grains. 

We initiated two delta-function velocity perturbations 
at t = 0 /is, one at each edge of the chain, both of magni¬ 
tude Vo = 9.899 X 10“^ mm//is, and directed into the 
chain. The resulting kinetic energy density plots are 
displayed in Figs, [^b) and id). The initial perturba¬ 
tions develop into SWs which travel toward the center 
of the chain, meeting precisely at the center due to the 
symmetry of the system. From Fig. [^b), it is evident 
that the central grain in the odd-numbered N = 39 
chain remains stationary during these collision processes, 
since the ISWs meet on the central grain itself. Never¬ 
theless, the large time scale associated with the rubber 
grains causes energy to remain in the vicinity of the soft 
grains for a noticeably long time. For the even-numbered 

= 40 grain system, the two ISWs meet at the interface 
between the central grains, thus there is no central cold- 
spot as there is in the three-soft grain system. This is 
clear from the kinetic energy density plot for the N = 40 
system (with four soft grains). Fig. ^d). This figure 
also illustrates that energy is localized near the four soft 
grains in this system for extended periods of time. The 
temporary energy localization can be attributed to the 


large time scale associated with the slow energy propa¬ 
gation in the soft impurities, and is also facilitated by 
the mutual repulsion between grains near the center of 
the chain. Also seen in Fig. [^b) is the steering of SSWs 
away from the central soft grains. Around 0.006 s, SSWs 
approaching the central grains are seen to turn away be¬ 
fore settling on grain 15 for a time. Likewise, at 0.009 s, 
two prominent SSWs stop on adjacent grains. This ki¬ 
netic energy localization is distinct from potential energy 
localization on the central grains, since being hard by 
comparison, these grains will transmit their energy very 
rapidly if they are in contact with another grain. How¬ 
ever, the arrival of a SW and a concomitant opening of a 
gap in the same direction of the SW prevents the imme¬ 
diate transmission of the energy and allows a temporary 
localization of kinetic energy at that grain location. 

To determine how breather-like these objects are, we 
analyze the frequency profiles via discrete cosine trans¬ 
forms (DCTs) of the square of the grain speed for one of 
the soft grains. These results are depicted in Fig. [To] for 
grain 19 in the three-soft grain system and for grain 21 
in the four-soft grain system. Here we have used the first 
100, 000 /is of the simulation. It is evident that in the 
three-soft grain system, a single frequency band develops 
centered between 4000 and 5000 Hz, whereas in the four- 
soft grain system, we see the development of multiple 
frequency bands: one wide band centered near 4000 Hz, 
and a more narrow band around 1000 Hz. The presence 
of more than one frequency band is a common charac¬ 
teristic of a breather m, suggesting that the four-soft 
grain system may be more breather-like of the two sys¬ 
tems considered here, and may therefore be particularly 
useful for locally trapping vibrational energy. Thus, our 
results show that this toy model may be one way to form 
breathers in granular chains. An extensive discussion on 
this aspect of the problem is under preparation. 


V. CONCLUSIONS 

We have performed here a detailed numerical analysis 
of material properties of granular chains, and the effects 
of these properties on the dynamics of, and journey to, 
the QEQ phase. A crucial finding of our study is that soft 
materials are characterized by an intrinsically large time 
scale associated with the propagation of energy, which 
causes a slow down of system dynamics in their vicinity. 
This slowing down can be exploited for multiple applica¬ 
tions. We see in our work the development of “sticky” 
walls in granular chains with soft boundaries, which are 
particularly useful for energy localization applications. 
We also presented a toy model that built on the results 
of the soft-wall systems, and seem to be promising for 
localizing energy at desired locations in a granular chain, 
with the resulting energy objects possessing properties 
similar to those of breathers. Beyond this, there is a 
natural beating response in the inhomogeneous soft-wall 
systems due to the interplay between large time scales of 
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FIG. 9. (Color online) A toy model for an energy focusing device. Here we have central soft grains composed of a soft material, 
such as rubber, embedded within a chain of stainless steel grains which are double the radius. To make the system as “soft” as 
possible, we include rubber walls at the edges of the stainless steel chain, (a) Three soft grains in an iV = 39 system, (c) four 
soft grains in an A" = 40 system. In each system, we initiate two delta-function velocity perturbations at t = 0 jus (one at each 
edge of the chain), directed toward the center of the chain, and then leave the systems to evolve in time. The corresponding 
kinetic energy density plots (with energy given on a logarithmic scale and normalized by the energy of the initial perturbation) 
are displayed, respectively, in (b) and (d). 
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FIG. 10. (Color online) Discrete cosine transforms of veloc¬ 
ities of soft grains in the toy model energy focusing devices. 
DCTs are performed over the first 100, 000 jus, illustrating the 
presence of a frequency band in the three-soft grain system 
with A = 39 grains, and multiple characteristic frequency 
bands in the four-soft grain system with A = 40 grains. Note 
that the low frequency components (/ < 200 Hz) have been 
removed from the DCT. 


the soft boundaries and small time scales of the harder 
grains. It is possible that this frequency response might 
be exploited by properly designed piezoelectric systems. 

One major finding of our work is that, since it takes 
longer for softer materials to rebound the grain-grain 
compression [35], incorporating softer materials at the 
boundaries of the granular chain causes a time delay 
in the rebounding of SWs, temporarily localizing en¬ 
ergy near the wall. This makes the soft walls appear 
“sticky” in dynamical simulations, and one noticeable 
consequence of this stickiness is the formation of larger 
SSWs. This was previously experimentally observed by 
Job et ah, who established that it is possible to increase 
the energy content of SSWs by softening the walls, at 
least for pB ^ 0.4 ^ 0.008 [19]. Moreover, since SSWs 
are produced slowly in the vicinity of the collision region, 
the soft-wall systems tend to have slightly hotter tem¬ 
peratures near a boundary for measurably long times, at 
least in dynamical simulations [29]. This was observed 
recently [35] via the formation of an energy gradient in 
the early time dynamics of an unloaded granular chain 
between asymmetric boundaries. 

We conclude by remarking on the equipartitioning of 
energy in the granular chain systems. In general, equipar¬ 
titioning of energy over time is facilitated by the nature 
of energy exchange in systems, which typically occurs be- 
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tween individual particles. If, on the other hand, inter¬ 
actions in a system do not permit the exchange of energy 
between individual particles, i.e. if the energy exchange 
must happen between small groups of particles m then 
each particle will not end up having more or less the 
same energy over time. The concept of temperature in 
such a system is not as well-defined, yet conceivably, the 
particles may still have a Gaussian distribution of veloc¬ 
ities (via Central Limit Theorem) and reach a state that 
is independent of initial conditions m- While this sys¬ 
tem would exhibit an equilibrium-like state, the absence 
of equipartitioning of energy - which manifests itself as 
large temperature fluctuations [32] - makes it distinct 
from the equilibrium state that is typically encountered 
in conventional solids, liquids, and gases. This is pre¬ 
cisely the nature of the QEQ state in the granular chain 
when subjected to a single delta function or short time 
perturbation at initiation. The QEQ state develops a 
sufficiently long time after an initial perturbation to the 


system. 
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Appendix: Material properties 

Table |T] lists the material parameters used in simula¬ 
tions, sorted from hard to soft according to pD where p is 
the material density. The prefactor between grains i and 
i + I is which is related to the Young’s modulus 

Y and Poisson’s ratio a of the grains by Equation 
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pD Young’s 
(mm//is)“^ 

modulus Poisson’s ratio 
(GPa) 

Density p Corresponding material 
(mg/mm^) 

0.00415 

1220 

0.2 

3.52 

Diamond 

0.00445 

720 

0.17 

2.2 

Quartz 

0.00805 

449.65 

0.207 

2.52 

Boron carbide 

0.0209 

249 

0.18 

3.58 

Magnesium oxide 

0.0479 

67 

0.2 

2.23 

Pyrex 

0.0551 

193 

0.305 

7.82 

Stainless steel 

0.0786 

11 

0.35 

0.657 

Red oak wood 

0.293 

97 

0.21 

19.816 

Plutonium 

0.367 

0.19 

0.5 

0.062 

Polyurethane-II 

0.450 

4.1 

0.39 

1.45 

PVC 

0.698 

0.1 

0.5 

0.062 

Poly ur et hane- III 

1.01 

0.069 

0.5 

0.062 

Polyurethane-I 

1.76 

1.46 

0.46 

2.17 

Teflon 

4.644 

0.55 

0.46 

2.16 Parylene teflon (PTFE) 

7.91 

0.032 

0.25 

0.18 

Cork 

26.5 

0.04 

0.4334 

0.87 

Polyolefin elastomer 

46.2 

0.03 

0.48 

1.2 

Rubber, hard 

52.4 

0.05 

0.49 

2.3 

Silicone rubber 


TABLE L Material properties used in simulations. 



